
/* Generated by Interface Builder */

#import "SignalProcessor.h"
#import <stdio.h>
#import <math.h>
#import <string.h>
#define SRATE 22050
#define PI 3.141592654782
#define SQRT_TWO 1.414213562
#define TWO_PI 6.283185309564

float sines[16384],cosines[16384];
int last_length = 0;

@implementation SignalProcessor

- window: (int) size array: (float *) array	//   call with defaults Hanning and size/2 centered
{
    return [self window: size array: array type: "Hanning" phase: 0];
}

- window: (int) size array: (float *) array type: (char *) type phase: (BOOL) phase
//	Window array of floats with window of type <"Triangle","Hanning","Hamming",
//		"Blackman3","Blackman4","Kaiser"> 
//		with window centered at zero (phase=TRUE) or size/2 (phase=FALSE)
{
    extern triangle(int size,float* array,BOOL phase),
    hanning(int size,float* array,BOOL phase),
    hamming(int size,float* array,BOOL phase),
    blackman3(int size,float* array,BOOL phase),
    blackman4(int size,float* array,BOOL phase),
    kaiser(int size,float* array,BOOL phase);

    if (!strcmp(type,"Triangle")) 		triangle(size,array,phase);
    else if (!strcmp(type,"Hanning")) 		hanning(size,array,phase);
    else if (!strcmp(type,"Hamming")) 		hamming(size,array,phase);
    else if (!strcmp(type,"Blackman3")) 	blackman3(size,array,phase);
    else if (!strcmp(type,"Blackman4")) 	blackman4(size,array,phase);
    else if (!strcmp(type,"Kaiser")) 		kaiser(size,array,phase);
// add your own windows here!!!
    else fprintf(stdout,"Windowing Rectangular\n");
    return self;
}

void triangle(int size,float* array,BOOL phase)
{
    int i;
    float w,win_frq;
    win_frq = TWO_PI / size;
    if (phase)
	for (i=0;i<size/2;i++)	{
	    w = size/2 - i;
	    array[i] *= w;
	    array[size-i] *= w;
	}
    else
	for (i=0;i<size/2;i++)	{
	    w = i;
	    array[i] *= w;
	    array[size-i-1] *= w;
	}
}

void hanning(int size,float* array,BOOL phase)
{
    int i;
    float w,win_frq;
    win_frq = TWO_PI / size;
    if (phase)
	for (i=0;i<size/2;i++)	{
	    w = 0.5 + 0.5 * cos(i * win_frq);
	    array[i] *= w;
	    array[size-i] *= w;
	}
    else
	for (i=0;i<size/2;i++)	{
	    w = 0.5 - 0.5 * cos(i * win_frq);
	    array[i] *= w;
	    array[size-i] *= w;
	}
}

void hamming(int size,float* array,BOOL phase)
{
    int i;
    float w,win_frq;
    win_frq = TWO_PI / size;
    if (phase)
	for (i=0;i<size/2;i++)	{
	    w = 0.54 + 0.46 * cos(i * win_frq);
	    array[i] *= w;
	    array[size-i] *= w;
	}
	else
	for (i=0;i<size/2;i++)	{
	    w = 0.54 - 0.46 * cos(i * win_frq);
	    array[i] *= w;
	    array[size-i] *= w;
	}
}

void blackman3(int size,float* array,BOOL phase)
{
    int i;
    float w,win_frq;
    win_frq = TWO_PI / size;
    if (phase)
	for (i=0;i<size/2;i++)	{
	    w = (0.42 + 0.5 * cos(i * win_frq) + 0.08 * cos(2 * i * win_frq));
	    array[i] *= w;
	    array[size-i] *= w;
	}
    else
	for (i=0;i<size/2;i++)	{
	    w = (0.42 - 0.5 * cos(i * win_frq) + 0.08 * cos(2 * i * win_frq));
	    array[size-i] *= w;
	    array[i] *= w;
	}
}

void blackman4(int size,float* array,BOOL phase)
{
    int i;
    float w,win_frq;
    win_frq = TWO_PI / size;
    if (phase)
	for (i=0;i<size/2;i++)	{
	    w = (0.35875 + 0.48829 * cos(i * win_frq) + 
			0.14128 * cos(2 * i * win_frq) + 0.01168 * cos(3 * i * win_frq));
	    array[i] *= w;
	    array[size-i] *= w;
	}
    else
	for (i=0;i<size/2;i++)	{
	    w = (0.35875 - 0.48829 * cos(i * win_frq) + 
			0.14128 * cos(2 * i * win_frq) - 0.01168 * cos(3 * i * win_frq));
	    array[i] *= w;
	    array[size-i] *= w;
	}
}

void kaiser(int size,float* array,BOOL phase)
{
    blackman4(size,array,phase);			//  I'm Cheating, see Harris paper, eq. 34
}

- logMag: (int) size array: (float *) f
{
    [self logMag: size array: f floor: -60.0];
    return self;
}

- log: (int) size array: (float *) f floor: (float) fl
{
    int i;
    float t1=0.0,t=0.0,t2;
    t2 = fl / 10.0;
    if (f[0] > t1)
        t1 = f[0];
    for (i=1;i<size;i++)	{
	if (f[i]>t1)
	   t1 = f[i];
    }
    for (i=0;i<size;i++)	{
    	t = log10(f[i]/t1);
	if (t<t2)
	   t = t2;
	f[i] = 1 - t/t2;
    }
    return self;
}

- log: (int) size array: (float *) f floor: (float) fl ceiling: (float) ceiling
{
    int i;
    float t1,t=0.0,t2;
    t2 = fl / 10.0;
    t1 = ceiling;
    for (i=0;i<size;i++)	{
    	t = log10(f[i]/t1);
	if (t<t2)
	   t = t2;
	f[i] = 1 - t/t2;
    }
    return self;
}

- logMag: (int) size array: (float *) f floor: (float) fl
{
    int i;
    float t1=0.0,t=0.0,t2=-12;
    t2 = fl / 10.0;
    f[0] = 2 * fabs(f[0]);
    if (f[0] > t1)
        t1 = f[0];
    for (i=1;i<size/2;i++)	{
    	t = f[i]*f[i] + f[size-i]*f[size-i];
	f[i] = t;
	if (t>t1)
	   t1 = t;
    }
    printf("%f\n",t1);
    for (i=0;i<size/2;i++)	{
    	t = log10(f[i]/t1);
	if (t<t2)
	   t = t2;
	f[i] = 1 - t/t2;
    }
    return self;
}

- logMag: (int) size array: (float *) f floor: (float) fl ceiling: (float) ceiling
{
    int i;
    float t1=0.0,t=0.0,t2=-12;
    t2 = fl / 10.0;
    f[0] = 2 * fabs(f[0]);
    for (i=1;i<size/2;i++)	{
    	t = f[i]*f[i] + f[size-i]*f[size-i];
	f[i] = t;
    }
    t1 = ceiling;
    for (i=0;i<size/2;i++)	{
    	t = log10(f[i]/t1);
	if (t<t2)
	   t = t2;
	f[i] = 1 - t/t2;
    }
    return self;
}

- magnitude: (int) size array: (float *) array magOut: (float *) mag
{
    int i;
    mag[0] = 2 * fabs(array[0]);
    for (i=1;i<size/2;i++)	{
        mag[i] = sqrt(array[i] * array[i] + array[size-i] * array[size-i]);
    }
    return self;
}

- square_mag: (int) size magnitudeIn: (float *) mag_in squareMagOut: (float *) square_mag_out
{
    int i;
    for (i=0;i<size;i++)  square_mag_out[i] = mag_in[i] * mag_in[i];
    return self;
}

- normalize: (int) size array: (float *) array;
{
    int i;
    float mx = 0.0;
    for (i=0;i<size;i++) if (fabs(array[i])>mx) mx=fabs(array[i]);
    if (mx!=0.0) for (i=0;i<size;i++) array[i] /= mx;
    return self;	
}

- dht: (float *) input_array output: (float *) output_array length: (int) length
{
    int i,j;
    float w;
    w = TWO_PI / length;
    for (i=0;i<length;i++)	{
        output_array[i] = 0.0;
	for (j=0;j<length;j++)	{
	    output_array[i] += input_array[j] * (cos(w * i * j) + sin(w * i * j));
	}
    }
    return self;
}

void make_sines(int length)
{
    int i;
    float freq,temp;
    if (length!=last_length)	{
        last_length = length;
        freq = 2.0 * PI / length;
        for (i=0;i<length;i++) {
            temp = freq * i;
	    sines[i] = sin(temp);
	    cosines[i] = cos(temp);
        }
    }
}

-  fhtRX2: (int) powerOfTwo array: (float *) array
{
    return self;		//   Not implemented yet
}

- fhtRX4: (int) powerOfFour array: (float *) array
{
    /*  In place Fast Hartley Transform of floating point data in array.
    	Size of data array must be power of four. Lots of sets of four 
	inline code statements, so it is verbose and repetitive, but fast. 
	A 1024 point FHT takes approximately 80 milliseconds on the NeXT computer
	(not in the DSP 56001, just in compiled C as shown here).

	The Fast Hartley Transform algorithm is patented, and is documented
	in the book "The Hartley Transform", by Ronald N. Bracewell.
	This routine was converted to C from a BASIC routine in the above book,
	that routine Copyright 1985, The Board of Trustees of Stanford University 	*/
	
    #define PI 3.141592654782
    #define SQRT_TWO 1.414213562
    #define TWO_PI 6.283185309564

    register int j=0,i=0,k=0,L=0;
    int n=0,n4=0,d1=0,d2=0,d3=0,d4=0,d5=1,d6=0,d7=0,d8=0,d9=0;
    int L1=0,L2=0,L3=0,L4=0,L5=0,L6=0,L7=0,L8=0;
    int n_over_d3;
    float r=0.0;
    int a1=0,a2=0,a3=0;
    float t=0.0,t1=0.0,t2 =0.0,t3=0.0,t4=0.0,t5=0.0,t6=0.0,t7=0.0,t8=0.0;
    float t9=0.0,t0=0.0;
    n = pow(4.0 , (double) powerOfFour);
    if (n!=last_length)	{
        make_sines(n);
	last_length = n;
    }
    n4 = n / 4;
    r = SQRT_TWO;
    j = 1;
    i = 0;
    while (i<n-1)	{
    	i++;
	if (i<j)	{
    		t = array[j-1];
		array[j-1] = array[i-1];
		array[i-1] = t;
    	}
    	k = n4;
    	while ((3*k)<j)	{
    		j -= 3 * k;
		k /= 4;
    	}
    	j += k;
    }
    for (i=0;i<n;i += 4) {
    	t5 = array[i];
	t6 = array[i+1];
    	t7 = array[i+2];
	t8 = array[i+3];
	t1 = t5 + t6;
    	t2 = t5 - t6;
    	t3 = t7 + t8;
    	t4 = t7 - t8;
	array[i] = t1 + t3;
	array[i+1] = t1 - t3;
	array[i+2] = t2 + t4;
	array[i+3] = t2 - t4;
    }
    for (L=2;L<=powerOfFour;L++)  {
	d1 = pow(2.0 , L+L-3.0);
	d2=d1+d1;
	d3=d2+d2;
	n_over_d3 = n / 2 / d3;
	d4=d2+d3;
	d5=d3+d3;
	for (j=0;j<n;j += d5)	  {
	     t5 = array[j];
	     t6 = array[j+d2];
	     t7 = array[j+d3];
	     t8 = array[j+d4];
	     t1 = t5+t6;
	     t2 = t5-t6;
	     t3 = t7+t8;
	     t4 = t7-t8;
	     array[j] = t1 + t3;
	     array[j+d2] = t1 - t3;
	     array[j+d3] = t2 + t4;
	     array[j+d4] = t2 - t4;
	     d6 = j+d1;
	     d7 = j+d1+d2;
	     d8 = j+d1+d3;
	     d9 = j+d1+d4;
	     t1 = array[d6];
	     t2 = array[d7] * r;
	     t3 = array[d8];
	     t4 = array[d9] * r;
	     array[d6] = t1 + t2 + t3;
	     array[d7] = t1 - t3 + t4;
	     array[d8] = t1 - t2 + t3;
	     array[d9] = t1 - t3 - t4;
	     for (k=1;k<d1;k++)	{
		  L1 = j + k;
		  L2 = L1 + d2;
		  L3 = L1 + d3;
		  L4 = L1 + d4;
		  L5 = j + d2 - k;
		  L6 = L5 + d2;
		  L7 = L5 + d3;
		  L8 = L5 + d4;
		  a1 = (int) (k * n_over_d3) % n;
		  a2 = (a1 + a1) % n;
		  a3 = (a1 + a2) % n;
		  t5 = array[L2] * cosines[a1] + array[L6] * sines[a1];
		  t6 = array[L3] * cosines[a2] + array[L7] * sines[a2];
		  t7 = array[L4] * cosines[a3] + array[L8] * sines[a3];
		  t8 = array[L6] * cosines[a1] - array[L2] * sines[a1];
		  t9 = array[L7] * cosines[a2] - array[L3] * sines[a2];
		  t0 = array[L8] * cosines[a3] - array[L4] * sines[a3];
		  t1 = array[L5] - t9;
		  t2 = array[L5] + t9;
		  t3 = - t8 - t0;
		  t4 = t5 - t7;
		  array[L5] = t1 + t4;
		  array[L6] = t2 + t3;
		  array[L7] = t1 - t4;
		  array[L8] = t2 - t3;
		  t1 = array[L1] + t6;
		  t2 = array[L1] - t6;
		  t3 = t8 - t0;
		  t4 = t5 + t7;
		  array[L1] = t1 + t4;
		  array[L2] = t2 + t3;
		  array[L3] = t1 - t4;
		  array[L4] = t2 - t3;
	     }
	}
    }		  
    return self;
}

@end
